function est_m = estimator_us_p4s3(est, estinfo)
% fixed effects
% s3: size cubic

nobs = size(est.YY,1);
ncity = length(est.xtval);
nyear = 6;

temp_y = est.YY;
temp_x = [est.XX, est.logDC, est.XX(:,14).^2, est.XX(:,14).^3];
% temp_yd = [year_dummy2, year_dummy3, year_dummy4, year_dummy5, year_dummy6];
temp_yd = est.YD;

islog = estinfo.islog;
if islog == 0
    temp_d = [est.DD, est.DD.^2, est.DD.^3, est.DD.^4];
elseif islog == 1
    temp_d = [log(est.DD),log(est.DD).^2, log(est.DD).^3, log(est.DD).^4];
elseif islog == 2
    temp_d = [log(1+est.DD), log(1+est.DD).^2, log(1+est.DD).^3, log(1+est.DD).^4];
end

xtflag = est.xtflag;
xtval = est.xtval;

xtflag2 = est.xtflag2;
xtval2 = est.xtval2;

xtflag3 = est.xtflag3;
xtval3 = est.xtval3;


% step 1: get rid of time effects
[d_y, m_y] = xt_demean(temp_y, [ones(nobs,1)], xtflag3, xtval3);
d_x     = xt_demean(temp_x, [ones(nobs,1)], xtflag3, xtval3);
d_d     = xt_demean(temp_d, [ones(nobs,1)], xtflag3, xtval3);

% step 2: get rid of heterogenous coeff on distance
dd_y = xt_demean(d_y, d_d, xtflag, xtval);
dd_x = xt_demean(d_x, d_d, xtflag, xtval);

% step 3: coefficients on common controls
dd_b = (dd_x'*dd_x)\(dd_x'*dd_y);

% step 4: coefficients on heterogenous controls
ddd_y = d_y-temp_x*dd_b;
ddd_x = d_d;
% ddd_b = (ddd_x'*ddd_x)\(ddd_x'*ddd_y);
[~, ddd_b] = xt_demean(ddd_y, ddd_x, xtflag, xtval);
[~,del_j] = xt_reg(ddd_b(:,1), ones(size(ddd_b,1),1), xtflag, xtval);
[~,del2_j] = xt_reg(ddd_b(:,2), ones(size(ddd_b,1),1), xtflag, xtval);
[~,del3_j] = xt_reg(ddd_b(:,3), ones(size(ddd_b,1),1), xtflag, xtval);
[~,del4_j] = xt_reg(ddd_b(:,4), ones(size(ddd_b,1),1), xtflag, xtval);

% step 5: alp
temp_resid = temp_y-temp_x*dd_b;
for i=1:1:size(temp_resid)
    temp_resid(i,:) = temp_resid(i,:) - temp_d(i,:)*ddd_b(i,:)';
end
[alp_jt_all,alp_jt] = xt_reg(temp_resid, ones(size(temp_resid,1),1), xtflag3, xtval3);
temp_sig2 = var(temp_resid-alp_jt_all);


% alp
alp_jt_mat  = nan*ones(length(xtval), length(xtval2));
alp_jt_mat2 = alp_jt_mat; %smooth version (extrapolate)
for j = 1:1:length(xtval)
    for t=1:1:length(xtval2)
        temp_xtval = xtval(j)+xtval2(t)/10000;
        temp_sel = xtval3 == temp_xtval;
        if sum(temp_sel) ~=0
            
            temp_alp = alp_jt(temp_sel,:);
            
            alp_jt_mat(j,t) = temp_alp;
            alp_jt_mat2(j,t) = temp_alp;
        else
            %alp_jt_mat2(j,t) = alp_j(j) + alp_t(j,t);
        end
    end
end

% store
est_m = [];
est_m.bet = dd_b(1:14);
est_m.bet2 = dd_b(16:17);

est_m.c10 = dd_b(15);

est_m.alp_j = alp_jt_mat2(:,1); %defined as the first year j
est_m.alp_jt_mat = alp_jt_mat;
est_m.alp_jt_mat2 = alp_jt_mat2;

est_m.del_j = del_j;
est_m.del2_j = del2_j;
est_m.del3_j = del3_j;
est_m.del4_j = del4_j;

est_m.model.MSE = temp_sig2;

% est_us2 = est_m;

% get the spatial function
est_m.sp = get_sp2(est_m, islog);


